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The advection and diffusion of a passive scalar is investigated for a map of the 2- 
torus. The map is chaotic, and the limit of almost-uniform stretching is considered. 
This allows an analytic understanding of the transition from a phase of constant 
scalar variance (for short times) to exponential decay (for long times). This transition 
is embodied in a short superexponential phase of decay. The asymptotic state in the 
exponential phase is an eigenfunction of the advection-diffusion operator, in which 
most of the scalar variance is concentrated at small scales, even though a large-scale 
mode sets the decay rate. The duration of the superexponential phase is proportional 
to the logarithm of the exponential decay rate; if the decay is slow enough then there 
is no superexponential phase at all. 
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A crucial problem involving fluids in the physical sciences is to understand 
the nature of mixing — its efficiency and thoroughness. Examples range from 
the mundane (cream in coffee), to the utilitarian (temperature in room), the 
industrial (mixing in chemical reactors), and the planetary (mixing of ozone 
in the extratropical stratosphere). If the flow is not turbulent, mixing can 
nevertheless be very efficient, due to a phenomenon called chaotic advection. 
In that case, the flow appears regular, but individual fluid trajectories are very 
complicated and lead to a stretching and folding action that greatly enhances 
mixing. Here we discuss mixing for a simple map, and show that a large-scale, 
coherent pattern is created that dominates the diffusive process. 

I. INTRODUCTION 

It has recently been suggested 1,2 that estimates of the decay rate of the variance of a 
passive scalar under the effect of advection and diffusion 3-6 do not yield satisfactory results 
when applied to some simple maps, such as the inhomogeneous baker's map. 7-9 This also 
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seems to be the case in laboratory experiments on periodic flows, 10 ' 11 where the decay rate 
is observed to be about an order of magnitude slower than the decay rate based on local 
arguments, such as the distribution of Lyapunov exponents. 12 Part of the reason for this 
is that in chaotic advection 13 (i.e., smooth flows with chaotic Lagrangian trajectories), far 
from the highly-turbulent regime, the presence of slowly-decaying eigenfunctions dominates 
the long-time decay rate. 1 ' 2 ' 14-16 (For the experiments, the presence of regular islands and 
barriers is also crucial, but we shall not address this complicated and poorly-understood 
issue here. It suffices to observe that the concentration field clearly attains an eigenfunction- 
like regime. 14 ) The existence of such eigenfunctions of the advection-diffusion operator was 
demonstrated convincingly via a numerical approach for the inhomogeneous baker's map. 1,2 
Sukhatme and Pierrehumbert 17 explained that the discrepancy is not due to a failure of the 
local approaches, but because they assume that the initial scale of variation of the passive 
scalar is much smaller than the system size. 

Here we propose to use a diffeomorphism of the 2-torus (an extension of Arnold's cat 
map 18 ) to further investigate aspects of the decay of variance and provide some analytical 
results. We find that, when the map is close to uniformly stretching, the decay rate is 
much faster than indicated by the distribution of Lyapunov exponents, as was also found in 
the inhomogeneous baker's map. 1 In Fereday et al. 1 and laboratory experiments, 12 a slower 
decay was also observed, but far from the uniformly-stretching (homogeneous) regime. 

The paper is organized as follows. In Section II we introduce the map and derive basic 
expressions for the effect of advection and diffusion on a passive scalar. We then analyze 
the superexponential (Section III) and exponential (Section IV) phases of diffusion. The 
spectrum of variance for the exponential eigenfunction is derived in Section V, followed by 
a discussion of the results in Section VI. 



II. ADVECTION DIFFUSION IN A MAP 

We consider a diffeomorphism of the 2-torus T 2 = [0, l] 2 , 

M(x)=M-x + <f>(x), (1) 

where M is a 2 x 2 nonsingular matrix with integer coefficients and (f>(x) is periodic in both 
directions with unit period. We choose M to have unit determinant, with an eigenvalue 
larger than one and the other less than one, so that even in absence of the <f> term M is still 
chaotic. Specifically, we take 

M=(?f); <f>(x) = - f^ 2 ™ 1 ) ; (2) 
V 1 V 2tt \sm2nx 1 J ' v ; 

so that M • x is the Arnold cat map and <p is a wave term usually associated with the 
standard map. The map M is area-preserving, and for K = the stretching of phase-space 
elements is uniform in space. The map is always chaotic (the largest Lyapunov exponent is 
positive). For smallK, there are no barriers to transport, such as islands, often encountered 
in realistic flows. 

We consider the effect of iterating the map and applying the heat operator to a scalar 
distribution 9^ l ~ 1 \x), 

6^(x) = 'K t {i - 1) (M-\x)), (3) 
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where e is the diffusivity, and the heat operator "K € and kernel h e are 

'K e e(x):=l h e (x- y)9(y)dy; h e (x) = V exp(27rifc • x - k 2 e). (4) 
J-P- ^ 



We Fourier expand 6^(x), 



so that (3) becomes 



0(0 (a.) = ^WeMfc-a, (5) 
^)( :r ) = ^T fc<J ^- 1 ), (6) 



with the transfer matrix, 

T kq := / exp (27ri (q ■ x - k ■ M(x)) - e q 2 ) dx. (7) 

Jo 2 

We may regard q as the "initial" wavenumber, and k as the "final" one, with a nonzero T kq 
denoting a transfer of concentration from q to k under one application of the map. 
For the form of the map given by (1) and (2), we have 

Tk q = e~ eq2 exp(2m(q- k-M) ■ x -[(k^ k 2 )K sm( Xl )) dx. (8) 
Jrp- 

The integral in x 2 gives a Kronecker delta, and the X\ integral is readily written as a Bessel 
function; we thus have 

T kq = e"^ 6 0tQ2 i Ql J Ql ((h + k 2 )K), Q := k ■ M - q, (9) 

where the Jq is a Bessel function of the first kind. 
In the absence of diffusion (e = 0), the variance 

a«:= f \eU(x)\ 2 dx = ^-> 4 l) H<f | 2 , (10) 

is preserved by (3) (we assume the spatial mean of 9 is zero), and for e > the variance 
decays (Fig. 1). We consider the case e <g 1, of greatest practical interest. For small K, 
there are three phases: (i) the variance is initially constant (if the initial scale of variation 
of the scalar concentration is well above the diffusive scale, as assumed here); (ii) it then 
undergoes a rapid superexponential decay; and (iii) it ultimately decays exponentially at 
a fixed rate, independent of e, as e — > 0. In the first phase, the map has not yet created 
gradients large enough for the small diffusion to act. In the second phase, there is a rapid 
exponential cascade to small scales and an associated exponential diffusion, leading to a 
superexponential decay. As the variance is depleted by diffusion, eventually the system 
settles into an eigenfunction that sets the exponential decay rate in the final phase. 

The existence of these three phases is well-known, 1_3 ' 19,20 but the exponential phase is 
the poorly understood, at least for the case of smooth flows and maps. We discuss the 
superexponential phase briefly in Section III, and in Section IV we describe the exponential 
phase. We will see that if the eigenfunction of the exponential phase decays slowly enough, 
then there is no superexponential phase at all. 
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FIG. 1: Decay of total variance for varying diffusivity e and K = 10~ 3 . The onset time of decay 
is logarithmic in the diffusivity, but the asymptotic exponential decay rate becomes independent 
of the diffusivity as e — > 0. The dashed curve shows the exact super exponential solution (K = 0) 
for e = 10~ 5 , and the dotted line is the single-mode value from Eq. (14). 

III. THE SUPEREXPONENTIAL PHASE 

Initially, the variance is essentially constant because the tiny diffusivity can be neglected. 
However, there is a cascade of the variance to larger wavenumbers under the action of M _1 
in (3). (In this phase, for small K, we can neglect the <f> term in (1), so that the map M 
is Arnold's cat map M • x.) This is the well-known "filamentation" effect in chaotic flows: 
the stretching and folding action of the flow causes rapid variation of the concentration 
across the folds. Thus, after a number of iterations %\ ~ 1 + (log e -1 / log A 2 ), 23 where A = 
(3 + \/5)/2 ~ 2.618 is the largest eigenvalue of M -1 , the diffusion can no longer be neglected. 
For e = 1CT 5 , we have i\ ~ 6 (this is always an overestimate). We now describe what happens 
to the variance after diffusion sets in. 

For small K and k, we have Jo ((fci + ki)K) ^> J\ ((ki + k,2)K), so we retain only 
the Q\ = term in the transfer matrix (9), 

T kq = e"^ 2 5 , Q + 0((h + k 2 ) 2 K 2 ) ■ (11) 

Hence, the nonvanishing matrix elements of T have k = q ■ M" 1 . If initially the variance 
is concentrated in a single wavenumber q (i.e., = unless k = q Q ), then after one 
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iteration it will all be in q • M -1 , after two in q ■ M -2 , etc. This amounts to the length 
of q being multiplied by a factor A > 1 at each iteration. But at each iteration the variance 
is multiplied by the diffusive decay factor exp(— 2eq 2 ), with q getting exponentially larger. 
The total variance is given by 

<r« = ,7 (0) exp(-2e ||q • M^ 1} || 2 ) ~ <7 (0) exp(-2e ||q || 2 A 2(i " 1} ), (12) 

so that the net decay is superexponential. The superexponential solution is represented by 
a dashed line in Fig. 1, with the solid line showing the numerical solution for the map M(a?). 
The superexponential solution is valid until about the ninth iteration. 24 We will revisit this 
breakdown of the solution in Section IV. 

It is to be noted that a more complicated initial condition also leads to superexponential 
decay, albeit with a less well-defined behavior because of the presence of several modes. Even 
an isotropic initial condition can be expected to have a superxponential phase: the averaging 
as performed in Antonsen, Jr. et al. 3 is problematic for a cat map in a periodic domain, 
because the slope of the stable (contracting) direction is irrational, and yet the wavevectors 
are confined to rational slopes. Hence, for finite k we cannot expect the averaging to hold. 
However, these difficulties are of a mathematical nature specific to the present problem and 
do not shed much light on a more general physical situation. 

IV. THE EXPONENTIAL PHASE 

In the superexponential phase we completely neglected the effect of the wave term in the 
map (1). We described the action as a perfect cascade to large wavenumbers, so that the 
variance was irrevocably moved to small scales and dissipated extremely rapidly. There can 
be no eigenfunction in such a situation, since the mode structure changes completely at each 
iteration. This direct cascade process dominates at first, but it is so efficient that eventually 
we must examine the effect of the the wave term, which is felt through the higher-order 
Bessel functions in the transfer matrix (9). 

Since the long-time exponential decay observed in the numerical results of Fig. 1 requires 
the existence of an eigenfunction, we may ask about the minimum requirement for this. 
Clearly if some scalar concentration is "left behind" in a given mode at each iteration, an 
eigenfunction will easily form. The question is then: Is it possible for the scalar concentration 
in a given wavenumber to be mapped back onto itself? This requires that the diagonal matrix 
element 

T qq = e-^U 0m i^J q2 {q 2 K) ) (13) 

be nonzero. We see from (13) that modes of the form q = (0 q 2 ) are mapped to themselves 
with a nonvanishing amplitude at each iteration: these are the modes that depend only on 
the X2 coordinate. This amplitude vanishes for K = 0, since q 2 ^ (the q = mode is 
preserved, and of no interest). 

For small K, the dominant Bessel function after J is J±, so the decay factor jj? for the 
variance is given by taking the magnitude of (13), 

A* = |T (0 i),(o i)| = e~ e Ji (K) = \K + 0{eK,K 2 ) . (14) 

Hence, for small K the decay rate is limited by the (0 1) mode. For e — * 0, the decay 
rate is independent of e. Figure 2 shows that the single-mode decay rate agrees very well 
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FIG. 2: Exponential decay rate log// 2 of the variance for e — > 0, as a function of K and. The 
triangles denote numerically calculated values, and the solid line is the small- -K" expression (14). 



with the numerical results even for K close to unity. In the inhomogeneous baker's map the 
nearly-superexponential limit is for a — > 1/2, where a is a parameter describing the inhomo- 
geneity of the map. For that case the transfer matrix scales in a manner analogous to here 
as a — > 1/2, but many more modes must be retained due to the presence of discontinuities: 
all the matrix coefficients decay as (1/2) — a, with none clearly dominating. The single-mode 
approximation is thus far less accurate. 

We can rule out the possibility that the decay is dominated by cycles that repeat after 
several iterations (that is, nonvanishing (T N ) qq for N > 1): such cycles must depend on 
higher-order Bessel functions that are small compared to J\(K). However, as K is made 
larger higher-order cycles become dominant and the situation is much more complicated. 

Now that the mechanism of exponential decay is understood (for small K) , we can go back 
and describe the condition for breakdown of the superexponential solution discussed at the 
end of Section III. The superexponential decay depletes the variance very rapidly until all 
that is left is variance in the exponentially decaying mode k := (0 1). The superexponential 
phase thus ends when the variance at large wavenumbers equals that in mode k . Assuming 
that the variance resides entirely in the fc mode initially, the condition for this is 

// 2 = exp(— e ||feo • M - ^ 2-1 ) || 2 ), (15) 

where \x is the decay factor of the variance in the fc mode, given by Eq. (14), and the right- 
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hand side is the superexponential solution (12). After substituting ||fco • M _ ^ 2_1 )|| ~ A* 2-1 , 
Eq. (15) must be solved numerically for i^. for K = 1(U 3 and e = 10~ 5 , we have %2 — 9.2. 
This is in fine agreement with the transition from superexponential to exponential in Fig. 1. 
If e <C 1, Eq. (15) has the approximate solution 

which gives %i ~ 8 for K = 10~ 3 , e = 10 -5 . Subtracting i± — 1 + log e -1 / log A 2 , the onset of 
the superexponential phase (Section III), we find that the duration of the superexponential 
phase is roughly 

log log /a" 1 

l 2-*i-— ; 7^—, (17) 

log A"* 

which is independent of e (at leading order), and has a weak dependence on the decay 
rate log/i. Unless \x is very small (recall that < ji < 1), the superexponential phase is very 
short. In fact, for log/i -1 < 1 the decay of the (0 1) mode is slow enough that there is 
no superexponential phase at all, as indicated by the negative right-hand side in (17). We 
can thus speculate that it is unlikely that the superexponential phase can be observed in 
experiments, since there \x tends to be close to unity. 

Note that e has to be extremely small for (17) to hold: for K = 10 -3 , e = 10 -5 , (17) 
gives %i — %\ ~ 1, whereas the unapproximated (numerical) result is %2 — %\ — 2.2. The error 
on (16) and (17) scales as log log e -1 . 



V. VARIANCE SPECTRUM OF THE EIGENFUNCTION 

The long-wavelength mode discussed in Section IV is the bottleneck that determines the 
decay rate (for small K). But this mode does not dominate the structure of the eigenfunction. 
In fact, a very small amount of the total variance actually resides in that bottleneck mode: 
the variance is concentrated at small scales. We now derive the variance spectrum of the 
eigenfunction. 

The variance is taken out of the (0 1) mode in the same manner as described in Sec- 
tion III: there is a cascade from that mode to larger wavenumbers through the action of M _1 . 
Neglecting the K term, the cascade proceeds from k = (0 1) to fei, fc 2 , &3, ... etc., as 

(0 1) -> (-1 2)->(-3 5)->(-8 13) (18) 

The k n become more and more aligned with the stable (contracting) direction of the map 
as we proceed down the cascade. The amplitude of the wavenumber is multiplied at each 
step by a factor A = (3 + v / 5)/2 ~ 2.618, the largest eigenvalue of Mr 1 . 

The exponential decay rate suggests that the scalar concentration is in an eigenfunction 
of the advection-diffusion operator. Assuming this to be the case, Fig. 3 illustrates the 
transfer of variance between modes for an iteration of the map. At each iteration, the 
eigenfunction property implies that the variance in each wavenumber is decreased by a 
uniform factor fi 2 < 1. This is illustrated by the vertical arrows in Fig. 3. The dashed 
arrows do not represent a direct transfer of variance, since for small K only the variance 
in the k mode is actually mapped back onto itself after one iteration (this is denoted by 
a solid vertical arrow). Rather, there is an effective (indirect) transfer occurring because of 
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FIG. 3: Schematic representation of the cascade of variance for an eigenfunction. The solid arrows 
represent a direct transfer of concentration, the dashed an "effective" transfer of amplitude \j? due 
to the eigenfunction property. In our aproximation, only the fco mode has a direct transfer of 
concentration to itself. 



the cascade (18): most of the variance in each mode is mapped to the next one down the 
cascade following the diagonal arrows in Fig. 3. 

The decrease in variance for each of the diagonal arrows is diffusive and is given by the 
factor v n = exp(— 2efc^). If we denote by crj^ := |#j?| 2 the variance in mode k n at the ith 
iteration, we have 

4l=f°V\ n = 0,l,..., (19a) 

4! = »n-i4:^ n = i,2,.... (i9b) 



These two recurrences can be combined to give 



v v m=0 / 



(20) 



where the relative variance in the nth mode is defined as a^(k n ) := cr^/a^}. The magnitude 
of the wavenumber is given by the exponential recursion, 

Ilfcj-Allfc^ll ||fc n || ~ A n ||fc || = A n , (21) 

which allows us to solve for n, 

n = log||fc n ||/logA (22) 

and rewrite (20) as 

a^\k n ) ~ ||fc n |r 21o ^ /losA exp {-2ekl/A 2 ) , (23) 

where we retained only the k\_ x term of the sum in (20) and used (21). The right-hand 
side of Eq. (23) for the relative variance does not (and should not if we really have an 
eigenfunction) depend on the iteration number, i, and depends only on n through k n . We 
thus let k n be a continuous variable k, with k = k (cos 9, sin 9), and drop i; from Eq. (23) 
we then have 

a(k, 9) = a(k) AT 1 5(9 - 9 ), (24) 
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FIG. 4: Spectrum function of relative variance after 12 iterations for K = 10 3 , e = 10 4 . The 
dashed line is the theoretical curve given by (25), and the triangles are numerical results. 



the spectrum of relative variance, with 

a(k) = k K exp (-2e k 2 /A 2 ) , C : = - log /V log A. (25) 

The factor k~ l 5(9 — 9 ) in (24) reflects the alignment of the vectors k with the stable 
(contracting) direction of M, which is at an angle 9 . We thus have essentially a one- 
dimensional spectrum. The spectrum function then satisfies 

a(k) = J a(k,9)kd9. (26) 

The spectrum function (25) is plotted in Fig. 4 and compared with numerical results 
for small K, showing excellent agreement. Since /i 2 < 1 and A > 1, we conclude that ( > 
always. This implies that there is more variance at the large wavenumbers than at the 
slowest-decaying mode k . The slope of the spectrum cr(k) is considerably shallower than 
the Batchelor k~ x spectrum, 21 consistent with the results for the baker's map. 1 This reflects 
the extreme efficiency of the cascade, a consequence of the nearly-uniform stretching, in that 
small scales are generated with great ease and the spectrum is therefore skewed towards large 
wavenumbers. 
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To know just how much more variance is at the large wavenumbers, we find the maximum 
of (25), which is at 

k m = A (C/2e) 1/2 , a(k m ) = k% e"< = k% // IogA . (27) 

The peak wavenumber thus scales as e -1 / 2 , the same scaling as the dissipation scale. 
From (27), the relative variance in that peak wavenumber scales as e _l ». The wavenum- 
ber k m gives an indication of the largest wavenumber that must be included in a numerical 
calculation to capture the decay of variance correctly. However, if the truncation size is 
smaller, the decay rate in the exponential phase is still captured properly, since it is deter- 
mined by the (0 1) mode. 

VI. DISCUSSION 

We summarize the three phases of chaotic mixing in smooth flows for the case of small 
diffusivity. In the first phase the variance is approximately conserved, and the chaotic 
flow (or map) creates large gradients in the scalar concentration through its stretching 
and folding action. This is usually called the stirring phase. In the second phase, the 
variance (that is, the squared-amplitude of each mode with the total mean subtracted) 
starts to decrease superexponentially, because the exponential cascade to small scales is 
compounded by the exponential efficiency of diffusion (Section III). This is the first of 
two mixing phases (superexponential and exponential), where diffusion plays an important 
role. This superexponential phase might not occur if the exponential decay rate of the 
slowest-decaying eigenfunction is slow enough. For very small diffusivity, the duration of 
the superexponential phase is independent of diffusivity. 

Unless the stretching is completely uniform, the superexponential phase comes to an end 
because though it rapidly depletes any variance contained in the small scales, some is left 
behind because of dispersion. What remains is the eigenfunction of the advect ion-diffusion 
operator with the largest eigenvalue (all eigenvalues have modulus less than one), which 
then decays exponentially. The decay rate of this eigenfunction is determined by its slowest- 
decaying part, in the present case the (0 1) mode (Section IV). The structure (spectrum) of 
this eigenmode is readily described as a balance between the eigenfunction property (modes 
are mapped to themselves with uniform amplitude) and a cascade to large wavenumbers 
(Section V). In the present case of a map with nearly uniform stretching, the spectrum of 
the eigenfunction has most of its variance concentrated at large wavenumbers, even though 
the small wavenumber mode (0 1) dictates the rate of decay. 

The decay rate of variance is outrageously fast in a map so close to being superexponential. 
Nevertheless, the manner in which the asymptotic regime is attained and the possibility of 
analytic results provide insight into the formation of the eigenfunction through the interplay 
of the slowest-decaying mode and the cascade to large wavenumbers. As K is made larger, 
the decay rates are more reasonable and a remnant of the mechanism presented here still 
applies. 

The decay rate in the present case is completely unrelated to the Lyapunov exponent or 
its distribution. For small K, the distribution of the Lyapunov exponent is peaked at log A 
and has a very narrow standard deviation. But here the asymptotic exponential decay 
rate is of order logK, so the decay becomes faster as K — > 0. This is due to the system 
being close to the uniform stretching (cat map) limit, which is unlikely to be the case in 



11 



physical situations. Any theory based on the distribution Lyapunov exponents cannot in 
this case predict the decay rate, since a global mode dominates. For the theory of Antonsen, 
Jr. et al. 3 , there is the further problem that, as in Fereday et al. 1 , averaging over angles 
is not possible here, since for small K the stable manifold and the gradient of the initial 
condition have a nearly constant angle with respect to each other. If the initial condition 
itself is taken as isotropic, then the irrationality of the slope of the contracting direction 
becomes problematic (Section III). 

Sukhatme and Pierrehumbert 17 point out that what determines the regime of decay (i.e., 
eigenfunction or local) is the scale of the initial scalar concentration. In our case, as we 
make that initial scale smaller we find the same asymptotic decay rate. This is due to a 
weak dispersion (due to the wave) from the large to the small wavenumbers which allows the 
system to develop its preferred (slowest-decaying) large-scale eigenfunction. The only way to 
get a faster rate is to completely remove the slowest-decaying eigenfunction from the initial 
condition, which never happens in practice. A smaller initial scale of variation does however 
lead to faster overall decay because its effect is to lengthen the initial superexponential scale. 
This is because the weak dispersion needs time to build the eigenfunction to an amplitude 
where it can rise above the other ambient modes. 

The large-scale eigenf unctions can lead not only to faster decay but also slower (as in Fere- 
day et al. 1 ), when compared to local, Lyapunov-exponent based approaches. 3 ' 4 ' 6 In both 
cases, it is the highly-ordered nature of the system (due to the large-scale, coherent nature 
of the initial scalar field and flow, but also to periodic boundary conditions and walls) that 
gives the discrepancy. We also observe a slower decay for larger K, but no analytical theory 
has yet been developed to adequately describe that regime. 

We observe numerically that as K is made large the spectrum of variance tends to concen- 
trate in small wavenumbers, possibly due to the presence of a strong dispersion competing 
with the direct cascade to small scales. 22 In that limit the cascade to large wavenumbers is 
no longer described by the linear part M of the map, so there is no clear separation between 
the eigenfunction property and the cascade. An investigation of the decay rate and spectrum 
in this large K, wave-dominated limit will be the subject of future work. 
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